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Abstract 

Migration of charged point defects triggered by the local random depolarization field is shown 
to plausibly explain aging of poled ferroelectric ceramics providing reasonable time and acceptor 
concentration dependences of the emerging internal bias field. The theory is based on the evaluation 
of the energy of the local depolarization field caused by mismatch of the polarizations of neighbor 
grains. The kinetics of charge migration assumes presence of mobile oxygen vacancies in the 
material due to the intentional or unintentional acceptor doping. Satisfactory agreement of the 
theory with experiment on the Fe-doped lead zirconate titanate is demonstrated. 

PACS numbers: 77.80.Dj,77.80.Fm,77.84.Dy,61.72.jd 
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I. INTRODUCTION 



Many ferroelectric materials exhibit gradual change of parameters with time under equi- 
librium external conditions, especially when acceptor doped. This phenomenon called ag- 
ing is known almost as long as ferroelectrics themselves^ but its mechanism is still highly 
disputed. Characteristic features of the degradation process are the decreasing dielectric 
constant and the fixed pattern of the polarization which hinders repolarization of the ma- 
terial. One of the first ideas of the aging mechanism was piling up of the space charge 
which pins the polarization configuration^^^^. Until recently, there was, however, no 
quantitative description of this mechanism which could allow comparison with experiment. 
An alternative and well elaborated concept of aging in acceptor doped ferroelectrics is the 
mechanism of defect dipole reorientation suggested by Arlt et al.->^ and supported in recent 
works^ 1 ^ 1 ^ 1 ^. This concept allowed reasonable explanation of the time and temperature 
dependencies of the most important parameters of aging, namely, of the emerging inter- 
nal bias field, E^,, and of the characteristic aging time, r. Nevertheless, this theory still 
seems to miss important features of aging concerning its dependence on the doping level. 
The orientation of the defect dipole due to the random walk of an oxygen vacancy about 
an acceptor defect is assumed to be a microscopic process independent of the other defect 
dipoles. Thus, the aging time appears to be independent of the doping level. The internal 
bias field proportional to the sum of independent contributions of the individual dipoles 
is expected in this theory to be proportional to the concentration of acceptor defects, Co- 
Experimentally, however, aging time is distinctly Cq dependent^^, and the internal bias field 
saturates with increasing concentration^^ well below Co — 1 mol%. For certain dopants the 
reason of this saturation could be the solubility limit in the host crystal as is stated to be 
the case for Fe-ions in PbZr x Ti!_ x 03 (PZT) ceramics^ 1 ^. The saturation of the bias field 
is observed, however, for virtually all acceptor dopants below Cq ~ 1 mol%^^ and seems to 
be a universal feature of acceptor doped ferroelectrics. 

As was recently shown, an alternative, charge defect migration mechanism can quanti- 
tatively explain essential features of aging in unpolarized ceramics^ 1 ^ 1 ^ as well as fatigue 
under a constant electric field^i. Following the latter concepts we advance in this work a 
model of aging in poled ferroelectric ceramics due to the depolarization-field driven charge 
migration. Theoretical results are then compared with model experiments on the PZT 
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ceramics with controlled Fe doping. 



II. MODEL OF FULLY POLARIZED FERROELECTRIC CERAMICS 

Ferroelectric ceramics are characterized by at least three sorts of randomness: random 
form of grains, their positions and random orientation of the crystal lattice inside the grains. 
We assume the last one to be the most important factor of randomness which can capture 
the main features of aging and fatigue in these systems, therefore only this kind of disorder 
will be considered in the following. 

We imagine the sample to consist of a regular cubic lattice of equal tightly contacted 
single-crystalline cubic grains of mesoscopic size R much larger than the lattice constant of 
the material. The grain edges are supposed to be aligned along the axes of the Cartesian 
coordinate system x, y, z as is shown in Fig. [TJ 

The sample is supposed to be sandwiched between plane electrodes located at z = ±L/2, 
where L 3> R, and polarized in a dc electric field substantially higher than the coercive 
field to the maximum possible spontaneous polarization in the z— direction. After that the 
voltage at the electrodes is set back to zero so that the remanent polarization P r in the 
z— direction remains. This state is considered as initial one for the process of aging at a 
fixed temperature and zero voltage maintained. 

Since the crystal structure of the grains is formed at temperatures far above the ferro- 
electric phase transition the crystal axes orientation in different grains is supposed to be 
absolutely arbitrary and independent. In the high electric field the polarization in every 
grain takes on the direction of one of the pseudocubic axes most close to the direction of the 
applied field (see in Fig. [1]). Note that the assumption that each grain has a homogeneous 
polarization is quite strong, as in reality grains always consist of a multitude of domains with 
different polarization directions. This is a well-known simplification that is often used in 
literature^ 1 ^. Thus, the vectors P s of local spontaneous polarization have the same magni- 
tude of P s and are randomly distributed within the cone defined by the polar angle 6 < 6 max 
with respect to the z— axis where 6 max = arcsin (y2/3) is the threshold angle introduced 
by Uchida et al.— . Further simplifying the problem we assume that the achieved maximum 
polarization in the field direction remains after setting the applied voltage back to zero. For 
calculation of thermodynamic characteristics one needs a procedure of statistical averaging 



3 



z 



■1 
111 




I 


11 

It 


J III 






1 1 1 




| - - - - - -1 


■n 


1 


111 



FIG. 1: Two-dimensional projection of polarization distribution in a fully polarized ferroelectric 
ceramics. Hatching pattern shows schematically the orientation of the crystal lattice. 



described in the next section. 



A. Configurational averaging 

Configurational averaging of local angle- dependent quantities over the ensemble of all 
possible random configurations, equal to the averaging over the sample volume, may be 
performed using the distribution function of a possible polarization direction in an arbitrary 
grain given by 

m<p) = %-, 0<9<n/4 and 

, vr/4 < 9 < 9 max (1) 

for < if < 2tt where (p and 9 are the azimuthal and polar angle in spherical coordinates 
associated with the above introduced cartesian coordinates and centered in the center of the 
chosen grain. The distribution function f(6,<p) is derived in Appendix lAl 

Let us introduce three-dimensional numeration (integer coordinates) of grains [n, k, m] 
associated with the cartesian coordinates so that the centers of the grains take positions 
(nR, kR,mR). Assuming independence of random angle variables in different grains the 
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distribution function for polarization directions in all grains reads 

(2) 

n,k,m 

where {6*j, ipi} denotes the manifold of spherical angles in all grains, while Q n ,k,m and (p n< k,m 
denote the angles in the grain with numbers [n, k, m]. 

Calculating an ensemble average of a quantity g({0i, <fi}): 

<9>=]jJ sin (W J &Pi9({0i,<Pi})F({0i,<Pi}) ( 3 ) 

one should take into account that all azimuthal and polar angle variables change in the same 
ranges as respective angles in Eq. ([I]). Thus, the ensemble average of the polarization along 
the z— axis, < P z >, amounts to 

P r = P s < cos {8 n k,m) >= P s ^—^ arcsin y/l/3, (4) 

7T 

and equals 0.831P S as in Refs.— & while the perpendicular polarization component averages 
to 

< P x >= P s < COS (9 n ,k,m) COS {(p n>k ,m) >= 0, (5) 

as well as < P y >= 0. The local polarization can be then conveniently decomposed in a sum 
of the mean and fluctuation polarizations as P s = P r + AP S where apparently < AP S >= 0. 



B. Mean values and variances of charges and fields 

Surface bound charge densities at the faces of a cubic grain with the number [n, k, m] 
located inside the bulk material are constant over the cubic faces and result from disconti- 
nuities of the respective normal components of the polarization AP S in the neighbor grains, 
namely, the charge density at the bottom face perpendicular to the axis z equals 

a n,k,m = P s (COS (0„,Jfc, m -l) - COS (9 nAm )) , (6) 

the charge density at the left face perpendicular to the axis x equals 

a n,k,m = P s [Sl n (#n-l,fc,m) COS (^ n -l,ifc,m) 

- Sin (0„,fc,m) COS ((p n ,k,rn)} , (7) 
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and the charge density at the left face perpendicular to the axis y equals 

°"n,fe,m = p s [sin (0 n ,fc-i,m) sin (<p n ,k-i,m) 

- sin (0njk, m ) sin ((p n ,k, m )] > ( 8 ) 

Configurational averaging of the above charge densities for internal grains with the function 
(j2J), equal to the averaging over the sample volume, results in vanishing mean values < 
(Tnt'm >= 0, which does not preclude the fact that local values ( l6fTf8]) are finite. For the top 
plane of the sample Eq. is not valid because there are no grains above the top grain layer. 
For that reason the nonzero mean value < a z n k m >= P r produced by the mean polarization 
P r prevails at the top plane of the sample, z = L/2. Similarly, < o z nkm >= —P r at the 
bottom plane of the sample, z = —L/2. 

Typical magnitudes of the charge densities at the internal grain faces ( 16TP71I81) are charac- 
terized by the variances of the respective charge densities and amount to 

° Z = \I< « fci J 2 > = 0.142- P, 
d*» = ^< KIJ 2 > = 0.547 • P s (9) 

The substantial difference between a z and a x ' y is explained by the fact that the direction 
of the polarization in the (x,y)— plane and, hence, variation of its azimuthal angle are 
not restricted while the polar angle is confined to the cone 6 < 6 max around the positive 
z— direction. 

Consistently with the splitting of polarization the total electric depolarization field can be 
decomposed as =< > +AE^ where the mean depolarization field < E^ > is created 
by the constant charge densities ±P r at the planes z = ±L/2, resulting from the mean 
polarization P r , and the fluctuation field AE^ is created by the fluctuation polarization AP S . 
In the considered poled state subject to aging, the field £ created inside the ferroelectric by 
the external sources (electrodes at z — ±L/2) exactly matches the mean depolarization field 
so that the mean total electric field < E >= £+ < E^ > vanishes together with the voltage 
at the electrodes, V — 0. This means that the mean surface charge densities at the top 
and bottom planes of the sample are precisely compensated by the opposite charge densities 
a = =FP r at the respective electrodes. 

Local depolarization field AE^ is connected by a nonlocal linear relation with the charge 
densities at all grain faces. To investigate statistic properties of this field we consider in the 
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following exemplary the field in the center of the central grain with numbers [0, 0, 0] which 
is representative for all internal grains of bulk ferroelectric ceramics. This field is given by 
the tensor relation 



(10) 
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J n,k,m p 

derived in the Appendix [0 Here indices a, (3 are introduced, taking on values x, y, z, and 
matrix elements T^ m are defined by the integrals over the cubic faces as is shown in the 
Appendix [B] The dielectric tensor of the ferroelectric medium is assumed to be of virtually 
cubic symmetry and is characterized by the relative permittivity e/. 

Typical magnitudes of the local electric field components are given by respective variances: 

/ p x 2 



V V T al3 ' <a p a 13 ' > 

/ j / j / j n,k,m n',k',m' n,k,m n' ,k' ,m' 

n,k,m n' ,k' ,m' f},/3' 

which are calculated in the Appendix [Bl The resulting typical field deviations are 



< (AE x d ' y ) 2 > = 3.523 



P, 



< (AE^ > = 0.898-^. (12) 
47reoe/ 

The preceding analysis demonstrates that local depolarization fields of the typical mag- 
nitude of P s /47reoe/ are present in bulk of the grains in the virgin state of poled ferroelectric 
ceramics. For the PZT materials with P s ~ 0.45 C/m 2 , and the lattice (high-field) value of 
the permittivity about e/ — 600 this field is about 2 ■ 10 7 V/m. 



C. Evolution of the electric field in the ceramics during aging 

If mobile charge carriers are present in the ceramics they have to be driven by the 
local electric field. Perovskite ferroelectrics are known to be semiconductors with pre- 
dominantly electronic or hole conductivity in reducing or, respectively, oxidizing atmo- 
sphere^£^>2£££i£L£L£iL£2.^ anc l a noticeable ionic contribution to the conductivity under the 
intermediate conditions. Specifically, PZT exhibits domination of ionic conductivity in a 
wide temperature range at atmospheric oxygen pressure^. In any case, the density of elec- 
tronic carriers in the samples quenched from high temperatures to the room temperature 
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is by many orders of the magnitude not sufficient to screen the surface bound charges in 
the ceramics 2 ^ 2 ^. In contrast to this, the density of oxygen vacancies, presumably equal to 
half the concentration of acceptor defects for electroneutrality reasons, is rather large since 
ferroelectrics are usually unintentionally acceptor doped^^ 1 ^ with small di- and trivalent 
cations which substitute for Zr 4+ or Ti 4+ . In this work the controlled Fe doping of PZT 
in the range of 0.1 — 1 mol% is maintained that results in the oxygen vacancies density of 
Co ~ 10 19 — 10 20 cm -3 , sufficient for screening of spontaneous polarization. An important 
question concerning the space-charge migration mechanism is: which part of the oxygen 
vacancies introduced by the acceptor doping is mobile. Measurements of ionic conductivity 
in Ca-doped barium titanate2§»2Z and in Ni-doped strontium titanate 2 ^ show a substantial 
increase of the conductivity with increasing doping which suggests that a significant part of 
the introduced vacancies is mobile. On the other hand, the ionic conductivity rises subpro- 
portionally with doping which may be the effect of defect association 2 ^. In the following we 
consider oxygen vacancies, the most mobile ionic species in perovskites, as a suitable agent 
for slow screening process. A possible effect of the subproportional increase of the mobile 
vacancy concentration with doping will be considered when adjusting theoretical results to 
the experimental data. 

During the migration of charge carriers they pile up at the charged grain faces to com- 
pensate the source of the electric field and thereby produce themselves the electric field Em- 
This process stops when the surface bound charges cr" fe m are outweighed by the emerging 
space charge. The charge redistribution inside the grains may be interpreted as a polariza- 
tion Pm, superimposed over the local polarization AP S . This is indeed a sort of migration 
polarization systematically studied experimentally by Okazaki in Refs.- 13 -^. In our model, 
however, this migration is induced by the fluctuation depolarization field AE^ which was 
not considered in the mentioned works by Okazaki. Typical thickness of space charge zones 
near the faces about h ~ a a /qco, with the oxygen vacancy charge equal twice the elemen- 
tary charge q and the vacancy concentration equal to half the acceptor concentration Co, 
amounts to 2 — 10 nm at cq = lmol%. This thickness is by two orders of magnitude less 
than the grain size R. This means that the space charge density <jm — hV~P m related to the 
migration polarization may be simply included in the total time dependent surface charge 
a nkm(t) = a rifcm(0) + °Af(£) where initial values at t — are given by the virgin charge 
densities (1617151) . 
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The influx of the screening charge to a charged face is driven by the local value of the 
electric field at the face which consists of the field generated by the face itself and that 
induced by the other faces. The latter contribution is continuous across the face and does 
not result in the net change of the surface charge. The former component is normal to the 
face at its both sides and equals <r" k m /2e e/ according to Gauss theorem. Thus the charge 
changes according to equation 

dt&n,k,m — ~ Ka n,k,m/ e O e f (13) 

where k = qfico is the conductivity due to oxygen vacancies with the mobility \l. Con- 
sequently, the complicated and random spatial distribution of the depolarization field in 
the sample remains unchanged during the charge migration while the magnitude of this 
field scales down coherently at all charged faces of the grains as c" fcm (0) exp (—t/r r ) where 
r r = e^Ef/n is the Maxwell- Wagner relaxation time. The latter seems to be the only char- 
acteristic time for aging mechanism through the charge migration assuming homogeneous 
background density of acceptor defects. 

In fact, there can be reasons for distribution of this time in a rather wide range. For 
example, if the activation energy of mobility is randomly distributed in some range this can 
result in quasilogarithmic time dependence of material parameters as in the case of aging 
in unpolarized ferroelectrics&2S. Another reason can be the field dependence of mobility as 
well as the complicated interplay of acceptor charge states and electronic state occupancies 
because of the band bending by the strong electric field near the charged faces of the grains 
as was suggested in photochemical studies of PZT-. The latter phenomena are still beyond 
the scope of our simple model. 

The relaxation time is determined by the concentration and mobility of the oxygen va- 
cancies. The former value is controlled by the acceptor concentration and can be reduced 
by the possible defect association as was mentioned before. The latter value is still highly 
disputed. Despite of the wide consensus concerning the migration energy barrier for oxygen 
vacancies about E a = 0.9 — 1.1 eV— i2L2&M the reported experimental values of the mobility 
at temperature T = 250° C range from 8 • 10~ 9 cm 2 /Vs , obtained by the simultaneous ther- 
moelectric power and conductivity measurements in Refs.— i2L2£ ; down to 2 • 1CT 14 cm 2 /Vs 
established from the diffusion depth profiles for oxygen tracers at the same temperature 
in Ref.— . Estimations at room temperature and the concentration of Co = 1 mol% result 
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accordingly in a wide range of possible aging times from 10 3 s to 4 • 10 8 s. In the Ref.— , how- 
ever, the strongly donor doped (1 — 4mol% of Nb 5+ ) PZT ceramics were studied in contrast 
to the acceptor doped ceramics of the Refs.— i2L2& which could have effect on the oxygen 
vacancy mobility. In the following we will rely on the ionic conductivity measurements of 
the acceptor doped PZT ceramics 26 . 



D. Thermodynamic analysis 

Considering the Gibbs free energy for ferroelectrics, Eq. (IC2I) of Appendix [Cj for the 
sample in the uniform external field exactly compensating the mean depolarization field, 
£ = — < Erf >, the contribution of the conductors ~ aV vanishes because of the zero 
voltage at the electrodes, V = 0, resulting in the initial energy before aging starts: 

G = F tt-- [ dVAP s AE d (14) 



2 _ 

where Q is the sample volume. The latter term here presents the contribution of the local 
fluctuations of the field and polarization to the total energy of the system and equals at 
t = 

AG + (0) = -Ifl J] < AP?AE% >= 0.649/2L (15) 

as is evaluated in the Appendix [Bl In the course of aging the local field and polarization are 
modified as AP S — > AP s + PAf(t) and AE^ — > AEd + Ejvf(£) and decrease as ~ exp (—t/r r ). 
Thereby the fluctuation contribution to the energy becomes 

AG+(t) = - l -J ^[AP s + P M (t)][AE d + E M (t)]. (16) 

This energy decreases with time ~ exp (— 2t/r r ) from the initial positive value, Eq. (fT5l) . 
driving the system to an energy minimum which means pinning of the given macroscopic 
polarization state. The strength of this pinning is characterized by the internal bias field 
which can be evaluated as follows. 

When reversing the external field to the opposite direction after the aging time t the 
local polarization and field change as P s — > —P s and E^ — > — E^. Thereby the fluctuation 
contribution to the energy changes to 

AGL(t) = ~f dV[~AP s + P M (t))[-AE d + E M {t)} (17) 
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where Pm(^) and Em(£) remain unchanged because the charge defects are too slow to follow 
the repolarization immediately. Taking into account that Pm(^) = — AP S [1 — exp {—t/r r )] 
and E M (t) = — AE d [l — exp {—t/r r )\ the difference in the energy of the opposite poled states 
amounts to 

AG_(*) - AG + (t) = J dV[AP s E M (t) + P M (t)AE d ] 

= 4AG+(0)[1 - exp (-t/r r )]. (18) 

Considering the increase in energy, which should be overcome by the repolarization, one can 
introduce an effective internal bias field as it was done by Arlt et al.— 

[AGL(t) - AG + (t)\ 



A 4^T |1_expH/Tr)1 (19) 



where A = 2.597 was calculated using Eq. ( Tl5j) . The maximum value of this field achieved 
at times t ^> r r is about 1.7 • 10 7 V/m which is in agreement with experimental estimations 
in Ref.- but exceeds by one order of the magnitude the values reported in Refs.— Very 
high theoretical value of the saturated bias field may follow from the basic hypothesis of the 
model on the absence of any correlations between polarizations in neighbor grains. That 
assumption allows appearance of unfavorable local configurations with high charge at the 
grain faces which can in reality be substantially depressed by local correlations. This fact 
will be considered below by adjusting the experimental curves. 

When applying to the aged sample the external field in the perpendicular direction as it 
is done in some experiments^ the fluctuation contribution becomes 

AG±(t) = -~ J dV[APj + P*(f)][AEtf + E M (t)\. (20) 

The fluctuation polarization in this state, AP^, is correlated neither with the initial fluctu- 
ation polarization AP S nor with the field AE^, therefore < AP^Eju(t) >= and, for the 
same reason, < P M (*)AE^ >= 0. On the other hand, < APjAEj >=< AP s AE d > since 
the system is macroscopically isotropic and this mean value should be direction independent. 
Consequently, the excess energy by the repolarization in perpendicular direction amounts to 

AG x (t) - AG+(t) = l -J dV[AV s E M {t) + V M {t)AE d ] 

= 2AG + (0)[l-exp(-t/r r )], (21) 
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and the internal bias field in the perpendicular direction is then equal half that in the initial 
polarization direction, E^(t) = (l/2)Ef b (t). 

We stress that the maximum magnitude of the bias field, Eq. (1191) . achieved at times 
t ^> r r , is defined by the energy of the fluctuation electrostatic field, AG + (0). This energy, 
Eq. (fT5|) . is proportional to the volume of the sample and, thus, is at least one order of 
the magnitude higher than the corresponding energy in the case of the unpolarized ferro- 
electrics^, which is accumulated near the charged domain faces. If the defect concentration 
in the bulk is so large that Co > c* = a a /qR and, hence, the thickness of the space-charge 
zone h <^ R, the whole energy of the fluctuation electrostatic field is virtually suppressed 
in the course of aging and, therefore, the maximum value of becomes concentration 
independent. In fact, the grain size R itself is doping dependent in the considered mate- 
rial which is typical for ferroelectric ceramics^ so that it is about 10 fim at cq — 0.1 mol% 
and decreases to 1 fim at Cq = 1 mol%. Thus, h <C R seems to be always the case. An 
apparent Co dependence of the measured may arise, however, due to the concentration 
dependence of the relaxation time r r . Thus, at the initial stage of aging, t <C r r , the bias 
field E ib (t) ~ t/r r ~ c . 

III. EXPERIMENT ON PZT 

To study the effect of doping on aging, PZT bulk samples with a zirconia to titania 
ratio of 54/46 were prepared with controlled contents of iron by a mixed oxide routed 
Undoped samples were used as well as samples doped with 0.1, 0.5 and 1 mol% of Fe. The 
powders were calcined, pressed into cylindrical bodies and sintered at 1050° C in air for 
6 hours. The bodies were cut into disc shaped samples with a wire saw (Well 2420, Well 
Diamantdrahtsaegen GmbH, Mannheim, Germany) and polished down to 3 /im grid size. 
After polishing the samples were relaxed in a furnace at a maximum temperature of 400° C 
for 2 hours to relieve mechanical stresses that were induced by the cutting and polishing 
process. Silver electrodes of 50 nm thickness were sputtered on the polished samples and a 
silver paste (Gwent Electronic Materials Ltd., Pontypool, U.K.,) was fired on afterwards at 
400° C. 

For the measurements it was necessary to start with fully unaged samples to ensure 
a statistical distribution of all defects and free charges. If the samples are cooled slowly 
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through the Curie temperature they age very fast and one can only obtain highly pinched 
hysteresis loops^. Therefore the samples were placed short circuited in a furnace at 450° C 
for 10 minutes before the first measurement. This leads to the necessary redistribution of 
the defects within the sample. After that they were quenched in pure water down to room 
temperature. By this the statistical charge configuration is "frozen" and the samples can 
easily be poled and measured in the unaged state. 

The polarization measurements were done in a Sawyer- Tower circuit with a 4.78 /xF test 
capacitance. The samples were placed between two soft springs in a silicon oil bath and 
the bipolar polarization hysteresis loops were measured at a frequency / = 55 mHz and a 
maximum field amplitude E = ±2kV/mm. The data were recorded by an oscilloscope. For 
each measurement step two bipolar cycles were applied. The first two cycles lead to the 
poling of the sample. After poling and measuring the samples were short circuited. Until 
the next measurement they were kept at room temperature in a sample box. Data points 
were taken after various time steps between 40 s and 10 6 s. 

The measured P-E loops were evaluated with respect to the internal bias field En,, to 
record the shift of the polarization curve along the field axis depending on the aging time^ 1 ^. 
This field was calculated as follows: 

E lb = -{E- +Et)/2 (22) 

where E~ is the coercive field when a negative external field is applied and is the coercive 
field on the positive field side. 

The internal bias field is presented as a function of time for different Fe contents in 
Fig. [21 Two common features of these curves are noticeable: they are virtually quasilinear 
on the logarithmic scale and all of them but one exhibit increasing slope. The first property 
is characteristic of aging and can be explained assuming wide distribution of relaxation 
times in the range r m j„ < r r < r max 9 , the second feature means then that the time span 
of measurements spreads well above r min but is still well below r max . Wide distribution of 
times r r may follow from a random distribution of activation energies Ej^ for which there 
are at least two physical reasons: various positions of the atomic cells with respect to the 
grain boundary, which can affect E a , and the mixture of regions with rhombohedral and 
tetragonal symmetry typical of the systems near the morphotropic phase boundary like the 
compound PbZro.54Tio.46O3 studied here. We assume in the following that the activation 
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FIG. 2: (Color online) Internal bias field as a function of time for PZT samples with different Fe 
contents. 

energy is a random variable with a Gaussian distribution g(E a ) ~ exp \—{E a — E a ) 2 /2s 2 } of 
the width s ~ 0.22 eV around the mean value E a ~ 1.12 eV. The internal bias field is then 
averaged over this distribution as 

poo 

E ib (t)= / dE a g(E a )E ib (E a ,t). (23) 
Jo 

where the bias field E ib (E a ,t) is given by Eq. (jl9]) with the relaxation time r r = e ef/qfic 
defined by the mobility 38 fi ~ T -1 exp (—E a /kT) with absolute temperature T and the 
Boltzmann constant k. 

By adjustment of the theoretical curves, Eq. (1231) . to the corresponding experimental 
ones for different doping concentrations two facts should be taken into account. Firstly, the 
above-mentioned grain size dependence on Co may have an effect on the maximum value of 
E ib . Indeed, it is known that larger grains in PZT often contain faults such as subgrains in 
their crystal structure^. These faults create regions with different polarization direction, but 
the correlation between these regions is higher than between neighboring independent small 
grains. This means higher energy of the fluctuation depolarization fields and, accordingly, 
higher E ib for smaller grains, i.e. for higher concentration cq. Secondly, accounting for the 
subproportional increase of the mobile vacancy concentration 28 with increasing Co and for the 
presumed solubility limit^^ for Fe in PZT at 0.5— 0.8mol% one should rather use an effective 
concentration of the mobile vacancies c e ff < Cq in the formula (|23|) . For those reasons, the 
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time (s) 

FIG. 3: (Color online) Fitting of the experimental time dependencies of the internal bias field of 
Fig. [2] (solid lines) with theoretical curves of Eq. (|23|) (dashed lines) is performed using the following 
field reduction factors and the effective concentrations of the mobile vacancies (see text): a = 0.141 
and c e fj = 0.15 mol% for c$ = 1.0 mol%, a = 0.127 and c e fj = 0.12 mol% for c$ = 0.5 mol%, 
a = 0.087 and c e jf = 0.09 mol% for cq = 0.1 mol%, and a = 0.045 and c e fj = 0.06 mol% for the 
nominally undoped PZT. 

reduction factor for the amplitude of the internal bias field, a < 1, and c e ff are used as 
fitting parameters when adjusting experimental results of Fig. [2] with Eq. (I2"5j) . As is seen 
in Fig. [3], satisfactory agreement is achieved by a ~ 0.05 — 0.1 and effective concentrations 
c eff ~ 0.15 mol% for c = 1.0 mol%, c e// ~ 0.12 mol% for c = 0.5 mol%, c e// ~ 0.09 mol% 
for Co = 0.1 mol% and c e // ~ 0.06 mol% for the nominally undoped material. Finally, the 
presented comparison of the charge migration model with experiments demonstrates that 
this mechanism may provide the observed magnitudes of the internal bias field and exhibit 
a reasonable agreement in the characteristic time and doping dependences. 

IV. CONCLUSIONS 

Even in the state with maximum spontaneous polarization of the ferroelectric ceramics 
the polarizations of individual grains cannot perfectly match each other all over the material 
because of random and independent orientation of the crystal lattice in different grains. That 
is why local bound charges at the grain surfaces and local depolarization field in the bulk are 
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present even if the mean depolarization field is exactly compensated in the electroded sample 
by zero voltage applied. We tried to capture this important feature of ferroelectric ceramics 
in the model of perfectly ordered cubic grains with random crystal lattice orientation in 
them. The frozen orientational disorder of the vectors of polarization in the grains can 
relax to the ordered state if mobile charges are present, which is the case in the acceptor 
doped (also unintentionally) ferroelectrics. Driven by the local depolarization field they 
can partly or completely compensate the fluctuation bound charge densities at the grain 
faces by building the thin space-charge zones near these faces. This is equivalent to effective 
alignment of the local polarizations in the direction of the total polarization accompanied by 
suppression of the fluctuation depolarization field. The described process leads system to the 
energy minimum where it resists its repolarization to another direction, i.e. exhibits aging 
of the poled ferroelectrics. Quantitative analysis of the relaxation kinetics in our model and 
its comparison to the model experiments with controlled acceptor contents show reasonable 
agreement in the magnitude of the internal bias field and its dependence on time and doping. 
By adjusting the experimental results with theoretical dependencies the activation energy 
of the vacancies was assumed smeared around the mean value of 1.1 eV, the concentration 
of the mobile vacancies c e // was taken reduced with respect to the acceptor concentration cq 
and the depression of the internal bias field by the possible correlation between the neighbor 
grains was taken into account. 

The suggested model assumes charge migration in the bulk of the grains which means, 
in fact, drift over the distance h, the thickness of the space charge zone. Another screening 
scenario may occur if the vacancies in the intergranular space are relevant to aging 5 . The 
vacancy diffusion along the grain boundaries is known to be two orders of the magnitude 
faster than the bulk diffusion 34 but the characteristic distance the charged defects should 
cover in this mechanism is given by the grain size R ^> h. The significance of the grain 
boundary contribution to screening seems therefore to increase for systems with smaller 
grains, i.e. by higher acceptor concentrations. This mechanism depends, however, on many 
still questionable factors, in the first instance, on the concentration of acceptors and oxygen 
vacancies in grain boundaries 1 ^ 1 ^ and their charge stated. 

Though the fluctuation depolarization field is present and causes charge migration in 
the whole volume of the polarized ceramics the charge defects pile up at the grain bound- 
aries making the considered effect a kind of surface phenomenon. This mechanism does not 
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preclude, however, other possible bulk mechanisms which can contribute to the aging of 
ferroelectric materials as, for example, the aforementioned defect dipole mechanism. For- 
mation of such Fe-0 and Cu-0 dipoles in the bulk lead titanate could be detected by means 
of electron paramagnetic resonance^ 2 -. It would be also interesting to prove by spectroscopic 
methods the presence of deficit or excess free vacancies and, respectively, Fe ions in space 
charge zones near the grain boundaries. This would be, however, not easy to detect because 
the volume share of these zones is about h/R ~ 10~ 3 — 10 -2 . 
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APPENDIX A: RANDOM DISTRIBUTION OF POLARIZATIONS IN GRAINS 

The distribution function f(9,f) can be derived following Uchida and Ikeda 23 . They 
considered a domain polarized in the positive z— direction and found a random distribution 
of the electric field, of the magnitude much higher than the coercive field, which is compatible 
with this direction of the polarization. Due to four-fold symmetry of the problem it is 
sufficient to consider the azimuthal region — 7r/4 < if < tt/4. A possible 90°-rotation of 
polarization from the z— to x— axis occurs when the projection of a unit field vector on 
the z— axis, cos#, becomes equal to the projection on the x— axis, cos 8' = sin 9 cos if. This 
determines the boundary of the area on the unit sphere 9(<p) = arctan(l/ cosy?). Thus the 
polar angle at this boundary changes between # = 7r/4at(/? = and 9 max = arctan (a/2) at 
if = tt/4. Within this area, the probability of different field directions which are compatible 
with the polarization along the z— axis is uniformly distributed with distribution function 
p(9,if) = 3/27r. Averaging of the field projection on the z— direction results in < cos# > = 
0.831, in accordance with Eq. (J3J). 

For the aims of this paper, we need a distribution function f(9, if) of the possible polar- 
ization orientations compatible with the strong electric field applied in the positive z— axis 
direction. In contrast to the function p(9, <f) this distribution is obviously not (^—dependent 
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since there is no special direction in the x — y plane. We assume here that the crystal lattice 
orientation is completely decoupled from the form and orientation of grains. The probability 
for the polarization to occur in the solid angle near the polar angle 9 is given by relation 

pw/4 

2irsin6d0f(6,(p) = 8 dtp sin 9 d9p(9,(p). (Al) 

Jip 

So long as the polar angle 9 is less than 7r/4 the integration over ip in the above integral goes 
from ifo = to 7r/4 which results in the constant value of f(0,ip). When tt/4 < < 9 max 
the integration goes from (po = arccos(cot#) to tt/4 which results in the function shown in 
Eq. p. 

APPENDIX B: CALCULATION OF THE LOCAL DEPOLARIZATION FIELD 
AND ITS MOMENTS 

The electric field ffTUl) is calculated by direct summation of contributions from all grain 
faces charged with densities (loTrTfgj) . Let us introduce notations n± = n± 1/2, k± = fc± 1/2, 
and m± = m ±1/2. Then the matrix elements T^ km in the formula (fTUj) are defined by the 
integrals over the cubic faces in dimensionless variables X = x/R,Y = y/R and Z = z/R: 

T^j\v£^ (nl+Y :- +Z2)y2 (B!) 
T »'L, = I dxj dZ {x2 + k2 _ + z2)3/2 (B2) 



n_ J m- 

+ r k + 



T n%, m ~ I dx I dY (x 2 + Y 2 + m 2 _) 3 / 2 



k- 



Y 



T ™> m= L dx L- dz {x 2 +k{+z 2 f/ 2 (B5) 

^ = /„_ dX L dY WTWT^- (B6) 

T nXm = I dX I dZ (X 2 + kl + Z 2 ) 3 / 2 ^ 



71- J 771 — 



m_ 



T^ m = I dX I dY {x2 + Y2 - +m2 _ )3/2 (B9) 



fe_ 
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After integration one obtains 



rpXX 



„./,.„i -= k +i n -) ~ 7(w_, k+, n-) 

- 7(777. + , n_) + 7(m_, k_,nJ), 



(BIO) 



where we introduced the function 

7(77-, k, m) = arctan 

and 



nk 



m\/n 2 + k 2 + m 2 



(Bll) 



T xy = In 

n,k,m 



+ In 



(m + + a/ 



77 2 + k 2 __ + TTlT 



(m + + \Jn 2 + J r k 2 _ + m 2 ^ 
+ a/t7+ + A; 2 + m 2 



m_ 



+ a/V + P + m 2 ^ 

For symmetry reasons the other components of the tensor T can be expressed as 



(B12) 



(B13) 



rpW _ rpXX 

n,k,m k,n,mi 

rpXZ __ rpXy 

n,k,m n,m,fc> 

rpyz _ rpxy 
n,k,m k,m,ni 

rpzy _ rpxy 

n,k,m m,k,n' 



rpzz _ rpXX 

n,k,m m,k,m 

rpyx _ rpxy 
n,k,m k,n,mi 

rpzx _ rpXy 

n,k,m m,n,fc> 



(B14) 
(B15) 
(B16) 
(B17) 



For calculation of the field component variances using Eq. (Hip the correlation functions 
of the charge densities have to be obtained first. They appear to be diagonal in Cartesian 
indices j3, (3' and involve only next neighbor indices n, k, m: 



< a n,k,m a n' ,k' ,m' > ~ Ps a 'i^k,k' ^m,m' 

< a n,k,m a n',k',m' > = P s a ^n,ri &m,m' 

x (25fc,fc' — Sk,k'-i — 8k,k'+i) 

X (2,5 mm > 3m,m' — 1 ^m.m'+l) 



(B18) 
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with constants 

, 1 2 
ai =< cos (6 nAm y >=- + — 7= = 0.701 

O 7TV3 

a 2 =< cos (9 n ,k,m) > 2 = s arcsin 2 (v/l/3) = 0.691 

a 3 =< sin [9 n k, m ) sin >= — - 

= 0.15. (B19) 

This is followed by expressions for the field component variances 

<(AE°f>= f-r^-Y (B20) 



4^606/^ 

E[~/7 r F ax f ryrriax rpax rpax \ 

L 3-*n,fe,m v n,k,m n—X,k,m n+l,k,m) 



X 

n,k,m 



_l_ „ T^ay (2T ay — T ay — T ay 1 

"i - u '3- i n,k,m \ n,k,m n,k—l,m n,k+l,m) 

4- (n — n \T az (r\rnaz _ rpaz _ rpaz 1 

l u l u 1) ± n,k,m \ Zj± n,k,m ± n,k,m-l ± n,k,m+l\ 

containing well converging series. The sums in the last formula are evaluated numerically 
resulting in formulas (|12p . 

The sum of correlation functions between the polarization and field arising in the Gibbs 
energy (|T5|) is position independent and can be evaluated, for example, for the central grain 
with n = k = m = using Eq. ([TO]) : 



AG + (0) = [a 3 (T- - T - ) (B21) 

o7reoe/ 



+°3 (^1,0 - ^0,0,o) + ~ a 2)(T 2 o i i - T 2 o i0 

resulting in Eq. ffToT) . 
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APPENDIX C: THERMODYNAMIC RELATIONS FOR A FERROELECTRIC 
SAMPLE IN AN EXTERNAL ELECTRIC FIELD 



An expression for the Gibbs free energy of a ferroelectric body subject to an external 
electric field is given by 36 

G = JdV\F - ^e E(e,,E) - EP S 

dV \F - ^ED - ^EP^ (CI) 

where E and D = P s + eoe r E are the total electrical field and the electrical displacement, 
respectively, P s denotes the local spontaneous polarization, and e r the relative permittivity 
tensor. The energy F depends on temperature and the material density only, while piezo- 
electric contributions are not considered. The integration in this expression is over the whole 
space excluding the volume of the conductors creating the external field. 

Using the potential presentation of the electric field E = — V0 and the fact that VD = 
due to the absence of the external charges one can transform ED = — V(0D) and then 
convert the integration of the respective term in Eqs. (IClj) to the integral over the surfaces 
of the conductors creating the external field. Since these surfaces are at constant poten- 
tials <pk arid the surface charge densities are given by the normal component of the electric 
displacement, D n , the Gibbs free energy can be transformed to 

G = ^-\H Q^k -\J dVEP,. (C2) 
k 

where Qk are the total charges at the conductors numerated by the index k, and JF depends 
only on temperature and the material density. In contrast to Eq. (1C1I) the integration in 
the last equation is over the volume of the dielectric body only since outside of the body 
the polarization vanishes. In equilibrium, the Gibbs free energy is minimum with respect to 
parameters involved 3 - when the potentials <fik at the conductors are kept constant. 
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